A cellular automaton model of gravitational clustering 
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Gravitational clustering of a random distribution of point masses is dominated by the effective 
short-range interactions due to large-scale isotropy. We introduce a one-dimensional cellular au- 
tomaton to reproduce this effect in the most schematic way: at each time particles move towards 
their nearest neighbours with whom they coalesce on collision. This model shows an extremely 
rich phenomenology with features of scale- invariant dynamics leading to a tree-like structure in 
space- time whose topological self-similarity are characterised with universal exponents. Our model 
suggests a simple interpretation of the non-analytic hierarchical clustering and can reproduce some 
of the self-similar features of gravitational N-body simulations. 
PACS numbers:89.75.Hc, 61.43.Hv, 98.80.-k 



The basic mechanism of large-scale pattern formation 
by gravity remains mainly elusive, in spite of elaborate 
studies in the past few decades [1]. Models of pattern 
formation, such as diffusion-limited aggregation (DLA) 
and its many variants [2], which have been remarkably 
successful in wide ranges of disciplines, do not appear to 
represent, even at a simplified level, gravitational clus- 
tering. A different class of irreversible aggregation mod- 
els, described by Smoluchowski coagulation equation and 
its various associated kernels, gives rise to a power-law 
distribution of the masses of the aggregates which are, 
however, randomly distributed in space [3-8]. 

Gravitational clustering, in general, leads to a power 
law distribution of the masses, known as the Press- 
Schechter mass function [9], but with a space distribu- 
tion which, at least up to some scale, is not homogeneous 
and can be described by fractal geometry [10-12]. The 
key difference between the aggregations one encounters 
in statistical physics and gravitational clustering is that 
in the former, forces are short-ranged, and the nonlin- 
ear dynamics, driven by collisions, erase the memory of 
the initial conditions rather fast. Gravity, on the other 
hand, is long-ranged with a deterministic evolution which 
strongly depends on the initial conditions even after the 
dynamics become nonlinear. 

Despite these notorious features of gravity, it has long 
been shown by Chandrasekhar [13] that for the Pois- 
son distribution of particle positions, gravity is effec- 
tively short-ranged : it can be well-described by a simple 
nearest-neighbours interactions. For the Poisson distri- 
bution of particles, the magnitude of the gravitational 
force, |F|, has a Holtsmark distribution [14] given by 
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where a = 4n(27rGm) 3 / 2 /15, which is a constant depend- 
ing on the mass of each particle, m, and the uniform 
number density n [13]. It was shown analytically and 
verified numerically that the 5/2 power-law tail of the 



Holtsmark distribution, is exactly due to the forces be- 
tween the nearest neighbours [13]. Even more remarkable 
is the fact that the distribution of the forces from the first 
neighbour approximation agrees with (1) over most of 
the range of |F| [13]. That the long-range gravitational 
force can be well-described by an effective short-range 
interaction may at first appear puzzling. However, the 
reason for this simplification is that the isotropy of Pois- 
son distribution is only broken at small scales where the 
granularity becomes important. 

As the system evolves, particle positions become cor- 
related, but the self-similar nature of gravitational clus- 
tering (see for example [9,15]) can ensure that at progres- 
sively larger scales the functional form of the Holtsmark 
distribution remains intact. Hence, under an appropriate 
renormalisation of the masses and the distances, deter- 
mined by the scale of granularity at a given time, gravi- 
tational interactions can be taken to be effectively short- 
ranged, even at times comparable to the dynamical time 
of the system. 

Inspired by these facts, we introduce an extremely 
simple one-dimensional automaton model which focuses 
precisely on this granular and self-similar nature of the 
gravitational clustering phenomenon. Unlike most previ- 
ous models such as DLA, Smoluchowski or 1-dimensional 
Burgers, our aggregation rule is position-dependent and 
independent of the masses and initial velocities of the 
aggregates. Unlike the conventional boolean cellular au- 
tomata, we distribute particles randomly on a line rather 
than on discrete lattice points. We focus mainly on the 
universal features of the distribution of the masses of the 
aggregates and on the space-time topological properties 
of the full aggregation history, i.e., on the universal char- 
acteristics of the "merger tree" of our model. 

The basic idea of this model, that the particles move 
towards their nearest neighbours, is implemented by the 
following algorithm. We distribute 10 5 particles ran- 
domly on a periodic line with a length of 10 5 units, so 
that the average density is equal to one. Particles move 
towards their nearest neighbours either by one unit at 
each time step or by half their separations, whichever 
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that is shorter. If two particles are closer than a lower 
threshold, which we set equal to 2 units, and, in addi- 
tion, are mutual nearest-neighbours they coalesce at the 
mid-point of their separation, conserving mass. This ag- 
gregation rule for a particle i at the position Xi (t) at time 
t can be written as, 



Xi(t + At) =Xi(t) +Min 
Xi(t + At) =Xi(t) +Max 
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If X r {t) < \ Xl {t)\ 
If X r (t) > \xi(t)\ 

(2) 



where x r (t) = x i+ i (t) — X{ (t) and x\ (t) = Xi-i(t) — X{ (t) . 
This process is illustrated in Fig. 1. 
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FIG. 1. A schematic illustration of the aggregation rule (2). 
The essence of our model is that particles move towards their 
nearest neighbours and coalesce on collision. 

As time progresses, masses and distances of the nearest 
neighbours rescale and the same rule, given by (2), is 
followed by more massive and farther apart aggregates. 
The continuation of this process will eventually lead to 
the total coagulation of the colloidal particles into one 
single mass. This aggregation mechanism traces a self- 
similar tree-like structure in space-time as shown in Fig. 
2. 

The tree structure of the aggregation process in space- 
time is a manifestation of topological self-similarity [16], 
which is a property of many branched structures such 
as river- networks [17] and bronchial trees [18] and can be 
quantified by various scaling exponents. A most common 
of these exponents is the Strahler index, s, given by 



(3) 



where P(r]) is the probability that a point in the network 
is connected to r] other points uphill, also known as the 
drainage area. The Strahler index, 8, is a measure of the 
bushiness of a branched structure and has an upper limit 
of 1.5 for river networks [17]. 

In our model, the function P(r?), given by (3), inter- 
prets as the probability that a newly-formed aggregate 
has a mass equal to r]. We observe a larger value for 
Strahler index (s = 2, given by the slope of the plot in the 
upper inset of Fig. 3) than is expected for river networks 
and which seems to be a characteristic of Cayley tree 



structures. A notable example of such a structure, also 
with exponent 2, has recently been observed for the prop- 
erty of internet connections [19]. The difference between 
our model and the conventional river-networks could be 
the stringent requirement of mass conservation here. 
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FIG. 2. The merger tree: the trajectories of 10 5 particles 
from a run of our simulation. The space-time merger tree 
exhibits the property of topological scaling, a concept which 
is relevant for branched structures and originally arose as a 
means of analysing river networks in two-dimensional space. 

Topological scaling, is believed to emerge from a self- 
similar growth process. Self-similar growth, or dynam- 
ical scaling, is a dominant feature of hierarchical grav- 
itational clustering [15], and is the fundamental reason 
for the emergence of distribution functions such as Press- 
Schechter mass function in cosmology [9] . An appropriate 
way to analyse dynamical scaling is, indeed, to study the 
mass distribution function n(m,t). Note that our func- 
tion differs from the usual number density by a constant 
factor which is given by the size of our system. 

It can be easily inferred from Fig. 2, and we have ver- 
ified numerically that the average mass, (m), and the 
mass variance, ((m 2 ) — (m) 2 ) 1 / 2 grow linearly with time. 
It is worth comparing this with the growth of average 
mass in one-dimensional Burgers, for example with uni- 
form initial velocities and positions, where an exponent 
of 2/3 has been obtained [20,21]. The linear growth of 
average mass also holds for Smoluchowski equation with 
a constant kernel [8] . The universality between our model 
and Smoluchowski arises inspite of the fact that particle 
trajectories are not Brownian here. In addition to this 
common feature, our model and Smoluchowski equation 
have similar asymptotic states. In both of these mod- 
els, clusters collide until all the mass falls into one final 
clump, whereas in 1-dimensional Burgers, the asymptotic 
state can contain many clumps with zero momentum. 

The distribution of the masses deviates slightly from 
a simple Gaussian for small masses where it develops a 
power-law as is shown in Fig. 4. The inset clearly demon- 
strates that a self-similar condensation process, rather 
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similar to the one observed in gravitational N-body sim- 
ulations (compare with Fig. 1 of [9]), has set in. In the 
process of coagulation, the shape of the mass spectrum 
seems to become fixed, and the curves move in parallel 
to the right (increasing aggregate mass). 



namely the power-law behaviour at small masses of the 
time-integrated mass distribution function, J n(m, t)dt ~ 
l/m l , is implied by the scaling relation (5), which fixes 
the value of the exponent I to I = 1. This value is also 
confirmed by our numerical simulations. 




FIG. 3. Topological scaling: the scaling, over many 
decades, of the probability distribution of the masses of the 
newly-borned aggregates, P(rj), (upper right inset) and of the 
corresponding cumulative probability, P> v , i-e. the probabil- 
ity that the mass of a newly-formed aggregate is larger than 
77, (main plot). The sketch in lower left inset illustrates the 
standard procedure of Horton-Strahler ordering and link am- 
plitude to find the Strahler index: P(rj) is the frequency of 
occurrence of a number shown in that sketch. 

The results presented so far provide strong evidence 
that the mass distribution function has the general scal- 
ing form: 



n(m,t) ~ m a t^ exp (— m 1 /t* 



(4) 



where the value of the exponents a,/?, 7 and A will be 
found in what follows. 

The conservation of the total mass in our model leads 
to the exponent identity 



/?+-(2 + a) 

7 



0. 



(5) 



In addition, the linear rate of growth of the average 
mass, (m) ~ t, leads to the further scaling identity, 
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(6) 



At this point one can already see the emergence of 
a topological scaling, namely a power-law behaviour in 
the time-integrated mass distribution function, at small 
masses. We comment that this is different from Strahler 
index we found in Fig. 3, since the latter refers only 
to newly-formed aggregates. This topological scaling, 



/ ♦ \ t=6 

/ V °- s 

/ \ 06 
/ 1 & 

f *f 5ft 0.4 

/ / * V 
i # / W 
/ / \\ t=7 


/ 

_ ' / / / 

/ / / / / / y"' s 


/ / / \ ' 6 26 44 62 60 98 116 134 
/ / / \ l*V t=9 rnDSS 

/ / □/ *s W \ \^?t=i2 

-ifc-n t T T i i *>-4^>-* a In -Tif-ii > * T T T ' 



FIG. 4. Dynamical scaling: self-similar growth of mass dis- 
tribution function n(m, t). Fits on the main plot are obtained 
with a simple Gaussian. The inset shows the self-similar 
growth of the fraction, /(m,t), of mass in objects of mass 
smaller than m, spanning over the significant part of the dy- 
namical time. 

A further test of our scaling identities (5) and (6), is 
provided by the evolution of the maxima, n max , of the 
mass distribution function, i.e., by the rate of decay of 
the peaks of the main plot in Fig. 4. Using the fact that 
the average mass grows linearly with time, in the scaling 
expression (4), we obtain n max (m,t) ~ l/£ 2 , which is 
again confirmed by our simulations. 

Thus, the scaling identities (5) and (6), reduce our 
scaling ansatz (4) to the self-similar form 
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The factor of 1/t 2 in (7) is a consequence of mass 
conservation and the linear growth rate of the aver- 
age mass, and has also been observed for the constant- 
kernel solution to Smoluchowski equation [8]: n{m,t) = 
At~ 2 exp(— 2m/t). This solution of Smoluchowski equa- 
tion, is obtained from our general solution (7) by the 
transformation t — > 2t and by using the following values 
of the exponents: a = 0, 7 = 1. We shall soon show that 
these exponents take different values in our model. 

To put our notation in accordance with that used in 
cosmology, we replace the time variable by the cut-off 
mass m*(t). As we have mentioned previously, our simu- 
lations show that a typical mass grows linearly with time, 
i.e. m* ~ t. The exponents a and 7 in expression (7) can 
be found by plotting n(m, t)m* 2 against the ratio m/m*, 
which we have done in Fig. 5. 
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FIG. 5. The mass function: the time-rescaled mass dis- 
tribution function, n(ra,£)ra* 2 , scatters around the fit (the 
continuous curve) for more than two decades in time. The 
characteristic cutoff mass, m*, grows linearly with time. 

We show in Fig. 5 that for over two decades in time 
scale the functional form of the mass function, given by 
the RHS of (7), is preserved. The fit in Fig. 5 sets the 
value of the unknown exponents in (7) to a ~ 3 and 7 ~ 2 
and we finally arrive at the approximate expression 

n(m,t) ^ ( e -(W^(t))» (g) 

V ^ m*(tf \m*(t)J V ; 

for the mass distribution function. We emphasis that, 
our solution is different from the constant-kernel solu- 
tion to Smoluchowski equation which is an special case of 
our general scaling solution (7), as previously noted. We 
comment that unlike some of the previous aggregation 
models used in cosmology [6,7], our distribution function 
cannot be formed from a white noise initial spectrum. It 
remains to be seen if, as for Smoluchowski with additive 
kernel, the introduction of a mass-dependent factor in our 
aggregation rule (2), would give rise to a Press-Schechter 
type mass function. 

We have also analysed the density-density correlation 
function which has a power-law behaviour with exponent 
— 1 at small scales, indicating that the mass is distributed 
on zero-dimensional objects, and a crossover to a con- 
stant value at large scales where it reminisces the initial 
Poisson distribution. The crossover length increases lin- 
early with time and the growth of the structures is dom- 
inated by the granular properties which are shifted from 
small to large scales. Thus, our model differs from the 
usual statistical models which generate fractals, where 
large-scale structures are built while substructures are 
preserved and not destroyed as is the case here. In this 
sense the present model does not generate asymptotic 
fractal structures in space. 

In conclusion, the seminal result of Chandrasekhar, 



that the long-range gravitational interactions between 
randomly distributed particles can be almost exactly re- 
placed by nearest-neighbour interactions, stimulated us 
to present a simple aggregation model which captures 
this profound feature of gravitating systems. We have 
shown that our model exhibits topological self-similarity 
over many decades of mass scale and dynamical scaling 
over many decades of temporal scale. These properties 
make it a simple and intuitive model for the study of 
gravitational hierarchical clustering. 
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